{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "4ed546c8-cbad-4721-ba6b-7d185222990e",
   "metadata": {},
   "source": [
    "# Example 4\n",
    "No source, Multiply-connected domain, Dirichlet/Neumann problem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "228a690d-545c-48b5-a2d5-6a3e40d7963b",
   "metadata": {},
   "outputs": [],
   "source": [
    "import calfem.geometry as cfg\n",
    "import calfem.mesh as cfm\n",
    "import calfem.vis as cfv\n",
    "import numpy as np\n",
    "from toolkits import gauss_line_int, gauss_tri_int\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9b47a234-100f-4115-a5ee-670e2b590fdd",
   "metadata": {},
   "source": [
    "## 1. Geometry"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "6e66b311-f596-421a-a94c-95d53734ac62",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Info    : GMSH -> Python-module\n",
      "235\n",
      "384\n",
      "384\n"
     ]
    }
   ],
   "source": [
    "# Number of holes\n",
    "nH = 2\n",
    "# Mesh size factor\n",
    "factor = 2\n",
    "\n",
    "# externel boundary ID\n",
    "dir_id_1 = 100\n",
    "dir_id_2 = 200\n",
    "dir_id_3 = 300\n",
    "neu_id_1 = 400\n",
    "# hole boundaryd ID\n",
    "hole_neu_id_1 = 500\n",
    "hole_neu_id_2 = 600\n",
    "hole_neu_id_3 = 700\n",
    "hole_neu_id_4 = 800\n",
    "\n",
    "# Geometry definition\n",
    "g = cfg.Geometry()\n",
    "\n",
    "r1 = 0.15\n",
    "r2 = r1/np.sqrt(2)\n",
    "points = [[0,1],[0.25,1],[0.625,1],[1,1],[1.375,1],[1.75,1],[2,1],  # 0-6\n",
    "          [0,0.75],[0.25,0.75],[0.625,0.75],[1,0.75],[1.375,0.75],[1.75,0.75],[2,0.75], # 7-13\n",
    "          [0,0.5],[0.25,0.5],[0.625-r1,0.5],[0.625-r2,0.5+r2],[0.625,0.5+r1],[0.625+r2,0.5+r2],[0.625+r1,0.5],[1,0.5], # 14-21\n",
    "          [1.375-r1,0.5],[1.375-r2,0.5+r2],[1.375,0.5+r1],[1.375+r2,0.5+r2],[1.375+r1,0.5],[1.75,0.5],[2,0.5], # 22-28\n",
    "          [0.625,0.5],[1.375,0.5], # 29-30 circle center\n",
    "          [0.625-r2,0.5-r2],[0.625,0.5-r1],[0.625+r2,0.5-r2],[1.375-r2,0.5-r2],[1.375,0.5-r1],[1.375+r2,0.5-r2], # 31-36\n",
    "          [0,0.25],[0.25,0.25],[0.625,0.25],[1,0.25],[1.375,0.25],[1.75,0.25],[2,0.25], # 37-43\n",
    "          [0,0],[0.25,0],[0.625,0],[1,0],[1.375,0],[1.75,0],[2,0]] # 44-50\n",
    "\n",
    "points = [[0,1],[0.4,1],[0.6,1],[0.8,1],[1.2,1],[1.4,1],[1.6,1],[2.0,1],  # 0-7\n",
    "          [0,0.7],[0.4,0.7],[0.6,0.7],[0.8,0.7],[1.2,0.7],[1.4,0.7],[1.6,0.7],[2.0,0.7], # 8-15\n",
    "          [0,0.5],[0.4,0.5],[0.6-r1,0.5],[0.6-r2,0.5+r2],[0.6,0.5+r1],[0.6+r2,0.5+r2],[0.6+r1,0.5],[0.8,0.5], # 16-23\n",
    "          [1.2,0.5],[1.4-r1,0.5],[1.4-r2,0.5+r2],[1.4,0.5+r1],[1.4+r2,0.5+r2],[1.4+r1,0.5],[1.6,0.5],[2,0.5], # 24-31\n",
    "          [0.6,0.5],[1.4,0.5], # 32-33 circle center\n",
    "          [0.6-r2,0.5-r2],[0.6,0.5-r1],[0.6+r2,0.5-r2],[1.4-r2,0.5-r2],[1.4,0.5-r1],[1.4+r2,0.5-r2], # 34-39\n",
    "          [0,0.3],[0.4,0.3],[0.6,0.3],[0.8,0.3],[1.2,0.3],[1.4,0.3],[1.6,0.3],[2,0.3], # 40-47\n",
    "          [0,0],[0.4,0],[0.6,0],[0.8,0],[1.2,0],[1.4,0],[1.6,0],[2,0]] # 48-55\n",
    "\n",
    "splines = [[0,1],[1,2],[2,3],[3,4],[4,5],[5,6],[6,7], # 0-6 top\n",
    "           [0,8],[1,9],[2,10],[3,11],[4,12],[5,13],[6,14],[7,15], # 7-14\n",
    "           [8,9],[9,10],[10,11],[11,12],[12,13],[13,14],[14,15], # 15-21\n",
    "           [8,16],[9,17],[9,19],[10,20],[11,21],[11,23],[12,24],[12,26],[13,27],[14,28],[14,30],[15,31], # 22-33\n",
    "           [16,17],[17,18],[22,23],[23,24],[24,25],[29,30],[30,31], #34-40\n",
    "           [16,40],[17,41],[34,41],[35,42],[36,43],[23,43],[24,44],[37,44],[38,45],[39,46],[30,46],[31,47], #41-52\n",
    "           [40,41],[41,42],[42,43],[43,44],[44,45],[45,46],[46,47], # 53-59\n",
    "           [40,48],[41,49],[42,50],[43,51],[44,52],[45,53],[46,54],[47,55],# 60-67\n",
    "           [48,49],[49,50],[50,51],[51,52],[52,53],[53,54],[54,55],# 68-74\n",
    "          ]\n",
    "long = [0,3,6,15,18,21,34,37,40,53,56,59,68,71,74]\n",
    "\n",
    "arcs = [[18,32,19],[19,32,20],[20,32,21],[21,32,22],\n",
    "        [22,32,36],[36,32,35],[35,32,34],[34,32,18], # 75-82\n",
    "        [25,33,26],[26,33,27],[27,33,28],[28,33,29],\n",
    "        [29,33,39],[39,33,38],[38,33,37],[37,33,25], # 83-90\n",
    "       ]\n",
    "\n",
    "surfs = [[0,8,15,7],[1,9,16,8],[2,10,17,9],[3,11,18,10],[4,12,19,11],[5,13,20,12],[6,14,21,13],\n",
    "         [15,23,34,22],[24,75,35,23],[16,25,76,24],[17,26,77,25],[26,27,36,78],\n",
    "         [18,28,37,27],\n",
    "         [29,83,38,28],[19,30,84,29],[20,31,85,30],[31,32,39,86],[21,33,40,32],\n",
    "         [34,42,53,41],[35,82,43,42],[81,44,54,43],[80,45,55,44],[36,46,45,79],\n",
    "         [37,47,56,46],\n",
    "         [38,90,48,47],[89,49,57,48],[88,50,58,49],[39,51,50,87],[40,52,59,51],\n",
    "         [53,61,68,60],[54,62,69,61],[55,63,70,62],[56,64,71,63],[57,65,72,64],[58,66,73,65],[59,67,74,66]]\n",
    "\n",
    "\n",
    "\n",
    "for x, y in points:\n",
    "    g.point([x, y])\n",
    "for i, s in enumerate(splines):\n",
    "    if i in long:\n",
    "        g.spline(s, el_on_curve=2*factor)\n",
    "    else:\n",
    "        g.spline(s, el_on_curve=factor)\n",
    "for c in arcs:\n",
    "    g.circle(c, el_on_curve=factor)\n",
    "\n",
    "for f in surfs:\n",
    "    g.struct_surf(f)\n",
    "\n",
    "bcs = {100:[7,22,41,60],\n",
    "       200:[68,69,70,71,72,73,74],\n",
    "       101:[0,1,2,3,4,5,6],\n",
    "       201:[14,33,52,67],\n",
    "       301:[75,76,77,78,79,80,81,82],\n",
    "       401:[83,84,85,86,87,88,89,90],\n",
    "       }\n",
    "\n",
    "for marker in bcs:\n",
    "    for i in bcs[marker]:\n",
    "        g.curve_marker(ID=i, marker=marker)\n",
    "cfv.drawGeometry(g)\n",
    "cfv.showAndWait()\n",
    "# print(len(points))\n",
    "\n",
    "mesh = cfm.GmshMesh(g, return_boundary_elements=True)\n",
    "mesh.elType = 2\n",
    "mesh.dofsPerNode = 1\n",
    "\n",
    "coords, edof, dofs, bdofs, elementmarkers, boundaryElements  = mesh.create()\n",
    "\n",
    "print(len(coords))\n",
    "print(len(edof))\n",
    "# print(elementmarkers)\n",
    "# print(boundaryElements)\n",
    "print(len(elementmarkers))\n",
    "\n",
    "# Display mesh\n",
    "cfv.draw_mesh(coords=coords, edof=edof, dofs_per_node=mesh.dofs_per_node,\n",
    "              el_type=mesh.el_type,filled=True,title=\"Mesh\")\n",
    "cfv.showAndWait()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "c194b6cf-1d47-4fd5-afd0-1eff4406605c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Info    : GMSH -> Python-module\n"
     ]
    }
   ],
   "source": [
    "el_on_curve = 5\n",
    "# Geometry definition\n",
    "g = cfg.Geometry()\n",
    "\n",
    "r1 = 0.15\n",
    "r2 = r1/np.sqrt(2)\n",
    "\n",
    "points = [[0,1],[0.4,1],[0.6,1],[0.8,1],[1.2,1],[1.4,1],[1.6,1],[2.0,1],  # 0-7\n",
    "          [0,0.7],[0.4,0.7],[0.6,0.7],[0.8,0.7],[1.2,0.7],[1.4,0.7],[1.6,0.7],[2.0,0.7], # 8-15\n",
    "          [0,0.5],[0.4,0.5],[0.6-r1,0.5],[0.6-r2,0.5+r2],[0.6,0.5+r1],[0.6+r2,0.5+r2],[0.6+r1,0.5],[0.8,0.5], # 16-23\n",
    "          [1.2,0.5],[1.4-r1,0.5],[1.4-r2,0.5+r2],[1.4,0.5+r1],[1.4+r2,0.5+r2],[1.4+r1,0.5],[1.6,0.5],[2,0.5], # 24-31\n",
    "          [0.6,0.5],[1.4,0.5], # 32-33 circle center\n",
    "          [0.6-r2,0.5-r2],[0.6,0.5-r1],[0.6+r2,0.5-r2],[1.4-r2,0.5-r2],[1.4,0.5-r1],[1.4+r2,0.5-r2], # 34-39\n",
    "          [0,0.3],[0.4,0.3],[0.6,0.3],[0.8,0.3],[1.2,0.3],[1.4,0.3],[1.6,0.3],[2,0.3], # 40-47\n",
    "          [0,0],[0.4,0],[0.6,0],[0.8,0],[1.2,0],[1.4,0],[1.6,0],[2,0]] # 48-55\n",
    "\n",
    "splines = [[0,1],[1,2],[2,3],[3,4],[4,5],[5,6],[6,7], # 0-6 top\n",
    "           [0,8],[1,9],[2,10],[3,11],[4,12],[5,13],[6,14],[7,15], # 7-14\n",
    "           [8,9],[9,10],[10,11],[11,12],[12,13],[13,14],[14,15], # 15-21\n",
    "           [8,16],[9,17],[9,19],[10,20],[11,21],[11,23],[12,24],[12,26],[13,27],[14,28],[14,30],[15,31], # 22-33\n",
    "           [16,17],[17,18],[22,23],[23,24],[24,25],[29,30],[30,31], #34-40\n",
    "           [16,40],[17,41],[34,41],[35,42],[36,43],[23,43],[24,44],[37,44],[38,45],[39,46],[30,46],[31,47], #41-52\n",
    "           [40,41],[41,42],[42,43],[43,44],[44,45],[45,46],[46,47], # 53-59\n",
    "           [40,48],[41,49],[42,50],[43,51],[44,52],[45,53],[46,54],[47,55],# 60-67\n",
    "           [48,49],[49,50],[50,51],[51,52],[52,53],[53,54],[54,55],# 68-74\n",
    "          ]\n",
    "\n",
    "longs = [0,3,6,15,18,21,34,37,40,53,56,59,68,71,74]\n",
    "\n",
    "arcs = [[18,32,19],[19,32,20],[20,32,21],[21,32,22],\n",
    "        [22,32,36],[36,32,35],[35,32,34],[34,32,18], # 75-82\n",
    "        [25,33,26],[26,33,27],[27,33,28],[28,33,29],\n",
    "        [29,33,39],[39,33,38],[38,33,37],[37,33,25], # 83-90\n",
    "       ]\n",
    "\n",
    "surfaces = [[0,8,15,7],[1,9,16,8],[2,10,17,9],[3,11,18,10],[4,12,19,11],[5,13,20,12],[6,14,21,13],\n",
    "         [15,23,34,22],[24,75,35,23],[16,25,76,24],[17,26,77,25],[26,27,36,78],\n",
    "         [18,28,37,27],\n",
    "         [29,83,38,28],[19,30,84,29],[20,31,85,30],[31,32,39,86],[21,33,40,32],\n",
    "         [34,42,53,41],[35,82,43,42],[81,44,54,43],[80,45,55,44],[36,46,45,79],\n",
    "         [37,47,56,46],\n",
    "         [38,90,48,47],[89,49,57,48],[88,50,58,49],[39,51,50,87],[40,52,59,51],\n",
    "         [53,61,68,60],[54,62,69,61],[55,63,70,62],[56,64,71,63],[57,65,72,64],[58,66,73,65],[59,67,74,66]]\n",
    "\n",
    "for x, y in points:\n",
    "    g.point([x, y])\n",
    "\n",
    "for i, s in enumerate(splines):\n",
    "    if i in longs:\n",
    "        g.spline(s, el_on_curve=2*el_on_curve)\n",
    "    else:\n",
    "        g.spline(s, el_on_curve=el_on_curve)\n",
    "\n",
    "for c in arcs:\n",
    "    g.circle(c, el_on_curve=el_on_curve)\n",
    "\n",
    "for f in surfaces:\n",
    "    g.struct_surf(f)\n",
    "\n",
    "\n",
    "bcs = {100:[7,22,41,60],\n",
    "       200:[68,69,70,71,72,73,74],\n",
    "       300:[75,76,77,78],\n",
    "       101:[0,1,2,3,4,5,6],\n",
    "       201:[14,33,52,67],\n",
    "       301:[79,80,81,82],\n",
    "       401:[83,84,85,86,87,88,89,90],\n",
    "       }\n",
    "\n",
    "for marker in bcs:\n",
    "    for i in bcs[marker]:\n",
    "        g.curve_marker(ID=i, marker=marker)\n",
    "        \n",
    "mesh = cfm.GmshMesh(g, return_boundary_elements=True)\n",
    "mesh.elType = 2\n",
    "mesh.dofsPerNode = 1\n",
    "\n",
    "coords, edof, dofs, bdofs, elementmarkers, boundaryElements  = mesh.create()\n",
    "cfv.draw_mesh(coords=coords, edof=edof, dofs_per_node=mesh.dofs_per_node,\n",
    "              el_type=mesh.el_type,filled=True,title=\"Mesh\")\n",
    "cfv.showAndWait()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "576437ef-7a8f-4f0a-81cc-8a2026a87c60",
   "metadata": {},
   "outputs": [],
   "source": [
    "def _cluster_edge(edges):\n",
    "    from collections import deque\n",
    "    clusters = []\n",
    "\n",
    "    e, _ = edges.popitem()\n",
    "    clusters.append(deque([e]))\n",
    "    \n",
    "    while len(edges) != 0:\n",
    "        quit = False\n",
    "        for i, j in edges:\n",
    "            for c in clusters:\n",
    "                if i == c[-1][-1]:\n",
    "                    c.append((i, j))\n",
    "                    edges.pop((i, j), None)\n",
    "                    quit = True\n",
    "                    break\n",
    "                if j == c[0][0]:\n",
    "                    c.appendleft((i, j))\n",
    "                    edges.pop((i, j), None)\n",
    "                    quit = True\n",
    "                    break\n",
    "            if quit:\n",
    "                break\n",
    "        else:\n",
    "            if len(edges) != 0:\n",
    "                e, _ = edges.popitem()\n",
    "                clusters.append(deque([e]))\n",
    "\n",
    "    return clusters\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "368bcf70-8ce8-4d49-898c-fda5ab3dcb9d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{(0, 54): 101, (54, 55): 101, (55, 56): 101, (56, 1): 101, (1, 57): 101, (57, 2): 101, (2, 58): 101, (58, 3): 101, (3, 59): 101, (59, 60): 101, (60, 61): 101, (61, 4): 101, (4, 62): 101, (62, 5): 101, (5, 63): 101, (63, 6): 101, (6, 64): 101, (64, 65): 101, (65, 66): 101, (66, 7): 101, (7, 74): 201, (74, 15): 201, (15, 99): 201, (99, 31): 201, (31, 124): 201, (124, 45): 201, (45, 145): 201, (145, 53): 201, (18, 159): 301, (159, 19): 301, (19, 160): 301, (160, 20): 301, (20, 161): 301, (161, 21): 301, (21, 162): 301, (162, 22): 301, (22, 163): 301, (163, 34): 301, (34, 164): 301, (164, 33): 301, (33, 165): 301, (165, 32): 301, (32, 166): 301, (166, 18): 301, (25, 167): 401, (167, 26): 401, (26, 168): 401, (168, 27): 401, (27, 169): 401, (169, 28): 401, (28, 170): 401, (170, 29): 401, (29, 171): 401, (171, 37): 401, (37, 172): 401, (172, 36): 401, (36, 173): 401, (173, 35): 401, (35, 174): 401, (174, 25): 401}\n"
     ]
    }
   ],
   "source": [
    "NeuMarker = [101,201,301,401]\n",
    "Neumann = {}\n",
    "for marker in NeuMarker:\n",
    "    for e in boundaryElements[marker]:\n",
    "        edge = tuple(np.array(e['node-number-list'])-1)\n",
    "        Neumann[edge] = marker    \n",
    "\n",
    "print(Neumann)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "a23d340d-c4f6-4a0e-a894-6ba3e989d08c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "59 1\n",
      "58 2\n",
      "57 3\n",
      "56 4\n",
      "55 5\n",
      "54 6\n",
      "53 7\n",
      "52 8\n",
      "51 9\n",
      "50 10\n",
      "49 11\n",
      "48 12\n",
      "47 13\n",
      "46 14\n",
      "45 15\n",
      "44 16\n",
      "43 17\n",
      "42 18\n",
      "41 19\n",
      "40 20\n",
      "39 21\n",
      "38 22\n",
      "37 23\n",
      "36 24\n",
      "35 25\n",
      "34 26\n",
      "33 27\n",
      "32 28\n",
      "31 29\n",
      "30 30\n",
      "29 31\n",
      "28 32\n",
      "27 33\n",
      "26 34\n",
      "25 35\n",
      "24 36\n",
      "23 37\n",
      "22 38\n",
      "21 39\n",
      "20 40\n",
      "19 41\n",
      "18 42\n",
      "17 43\n",
      "16 44\n",
      "15 45\n",
      "14 46\n",
      "13 47\n",
      "12 48\n",
      "11 49\n",
      "10 50\n",
      "9 51\n",
      "8 52\n",
      "7 53\n",
      "6 54\n",
      "5 55\n",
      "4 56\n",
      "3 57\n",
      "2 58\n",
      "1 59\n",
      "[deque([(174, 25), (25, 167), (167, 26), (26, 168), (168, 27), (27, 169), (169, 28), (28, 170), (170, 29), (29, 171), (171, 37), (37, 172), (172, 36), (36, 173), (173, 35), (35, 174)]), deque([(166, 18), (18, 159), (159, 19), (19, 160), (160, 20), (20, 161), (161, 21), (21, 162), (162, 22), (22, 163), (163, 34), (34, 164), (164, 33), (33, 165), (165, 32), (32, 166)]), deque([(0, 54), (54, 55), (55, 56), (56, 1), (1, 57), (57, 2), (2, 58), (58, 3), (3, 59), (59, 60), (60, 61), (61, 4), (4, 62), (62, 5), (5, 63), (63, 6), (6, 64), (64, 65), (65, 66), (66, 7), (7, 74), (74, 15), (15, 99), (99, 31), (31, 124), (124, 45), (45, 145), (145, 53)])]\n"
     ]
    }
   ],
   "source": [
    "print(_cluster_edge(Neumann))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "10d1592c-2e93-4e35-a559-d4082314fd69",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{}\n"
     ]
    }
   ],
   "source": [
    "print(Neumann)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "b59c5b9f-5e61-4711-82a0-d4f6f1ab0af5",
   "metadata": {},
   "outputs": [],
   "source": [
    "g.point([0.0, 0.0], ID=0)\n",
    "g.point([1.0, 0.0], ID=1)\n",
    "g.point([1.0, 1.0], ID=2)\n",
    "g.point([0.0, 1.0], ID=3)\n",
    "g.point([0.4, 0.4], ID=4)\n",
    "g.point([0.6, 0.4], ID=5)\n",
    "g.point([0.6, 0.6], ID=6)\n",
    "g.point([0.4, 0.6], ID=7)\n",
    "g.spline([0, 1], el_on_curve=factor, marker=dir_id_1, ID=0)\n",
    "g.spline([1, 2], el_on_curve=factor, marker=neu_id_1, ID=1) \n",
    "g.spline([2, 3], el_on_curve=factor, marker=dir_id_2, ID=2) \n",
    "g.spline([3, 0], el_on_curve=factor, marker=dir_id_3, ID=3)\n",
    "g.spline([4, 5], el_on_curve=factor/4, marker=hole_neu_id_1, ID=4)\n",
    "g.spline([5, 6], el_on_curve=factor/4, marker=hole_neu_id_2, ID=5)\n",
    "g.spline([6, 7], el_on_curve=factor/4, marker=hole_neu_id_3, ID=6)\n",
    "g.spline([7, 4], el_on_curve=factor/4, marker=hole_neu_id_4, ID=7)\n",
    "g.surface([0, 1, 2, 3],[[4,5,6,7]])\n",
    "cfv.drawGeometry(g)\n",
    "cfv.showAndWait()\n",
    "\n",
    "# g.point([0.0, 0.0], ID=0)\n",
    "# g.point([1.0, 0.0], ID=1)\n",
    "# g.point([1.0, 1.0], ID=2)\n",
    "# g.point([0.0, 1.0], ID=3)\n",
    "# g.point([0.5, 0.5], ID=4)\n",
    "# g.point([0.6, 0.5], ID=5)\n",
    "# g.point([0.5, 0.6], ID=6)\n",
    "# g.point([0.4, 0.5], ID=7)\n",
    "# g.point([0.5, 0.4], ID=8)\n",
    "# g.spline([0, 1], el_on_curve=factor, marker=dir_id_1, ID=0)\n",
    "# g.spline([1, 2], el_on_curve=factor, marker=neu_id_1, ID=1) \n",
    "# g.spline([2, 3], el_on_curve=factor, marker=dir_id_2, ID=2) \n",
    "# g.spline([3, 0], el_on_curve=factor, marker=dir_id_3, ID=3)\n",
    "# g.circle([5, 4, 6], el_on_curve=factor/4, ID=4)\n",
    "# g.circle([6, 4, 7], el_on_curve=factor/4, ID=5)\n",
    "# g.circle([7, 4, 8], el_on_curve=factor/4, ID=6)\n",
    "# g.circle([8, 4, 5], el_on_curve=factor/4, ID=7)\n",
    "# g.surface([0, 1, 2, 3],[[4,5,6,7]])\n",
    "# # g.structuredSurface([0, 1, 2, 3])\n",
    "# Display mesh\n",
    "\n",
    "\n",
    "## 2. Meshing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "07931634-1ddc-4c54-9b23-e70b883c5d48",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Info    : GMSH -> Python-module\n",
      "[[  4 247]\n",
      " [247 246]\n",
      " [246 245]\n",
      " [245 244]\n",
      " [244 243]\n",
      " [243 242]\n",
      " [242 241]\n",
      " [241 240]\n",
      " [240 239]\n",
      " [239 238]\n",
      " [238 237]\n",
      " [237   7]\n",
      " [  7 236]\n",
      " [236 235]\n",
      " [235 234]\n",
      " [234 233]\n",
      " [233 232]\n",
      " [232 231]\n",
      " [231 230]\n",
      " [230 229]\n",
      " [229 228]\n",
      " [228 227]\n",
      " [227 226]\n",
      " [226   6]\n",
      " [  6 225]\n",
      " [225 224]\n",
      " [224 223]\n",
      " [223 222]\n",
      " [222 221]\n",
      " [221 220]\n",
      " [220 219]\n",
      " [219 218]\n",
      " [218 217]\n",
      " [217 216]\n",
      " [216 215]\n",
      " [215   5]\n",
      " [  5 214]\n",
      " [214 213]\n",
      " [213 212]\n",
      " [212 211]\n",
      " [211 210]\n",
      " [210 209]\n",
      " [209 208]\n",
      " [208 207]\n",
      " [207 206]\n",
      " [206 205]\n",
      " [205 204]\n",
      " [204   4]]\n",
      "[[  5 214]\n",
      " [214 213]\n",
      " [213 212]\n",
      " [212 211]\n",
      " [211 210]\n",
      " [210 209]\n",
      " [209 208]\n",
      " [208 207]\n",
      " [207 206]\n",
      " [206 205]\n",
      " [205 204]\n",
      " [204   4]]\n",
      "[[  6 225]\n",
      " [225 224]\n",
      " [224 223]\n",
      " [223 222]\n",
      " [222 221]\n",
      " [221 220]\n",
      " [220 219]\n",
      " [219 218]\n",
      " [218 217]\n",
      " [217 216]\n",
      " [216 215]\n",
      " [215   5]]\n",
      "[[  7 236]\n",
      " [236 235]\n",
      " [235 234]\n",
      " [234 233]\n",
      " [233 232]\n",
      " [232 231]\n",
      " [231 230]\n",
      " [230 229]\n",
      " [229 228]\n",
      " [228 227]\n",
      " [227 226]\n",
      " [226   6]]\n",
      "[[  4 247]\n",
      " [247 246]\n",
      " [246 245]\n",
      " [245 244]\n",
      " [244 243]\n",
      " [243 242]\n",
      " [242 241]\n",
      " [241 240]\n",
      " [240 239]\n",
      " [239 238]\n",
      " [238 237]\n",
      " [237   7]]\n",
      "Dirichlet !\n"
     ]
    }
   ],
   "source": [
    "mesh = cfm.GmshMesh(g,return_boundary_elements=True)\n",
    "mesh.elType = 2          # Degrees of freedom per node.\n",
    "mesh.dofsPerNode = 1     # Factor that changes element sizes.\n",
    "# mesh.elSizeFactor = 0.5 # Element size Factor\n",
    "\n",
    "coords, edof, dofs, bdofs, elementmarkers, boundaryElements  = mesh.create()\n",
    "\n",
    "#############################\n",
    "#                           #\n",
    "# Elements, Edges and Nodes #\n",
    "#                           #\n",
    "#############################\n",
    "element = edof - 1\n",
    "node = coords\n",
    "edge = dict()\n",
    "for i, e in enumerate(element):\n",
    "    edge[(e[0], e[1])] = i\n",
    "    edge[(e[1], e[2])] = i\n",
    "    edge[(e[2], e[0])] = i\n",
    "\n",
    "######################\n",
    "#                    #\n",
    "#   Externel edges   #\n",
    "#                    #\n",
    "######################\n",
    "dir_edge_1 = np.array([edge['node-number-list'] for edge in boundaryElements[dir_id_1]])-1\n",
    "dir_edge_2 = np.array([edge['node-number-list'] for edge in boundaryElements[dir_id_2]])-1\n",
    "dir_edge_3 = np.array([edge['node-number-list'] for edge in boundaryElements[dir_id_3]])-1\n",
    "\n",
    "neu_edge_1 = np.array([edge['node-number-list'] for edge in boundaryElements[neu_id_1]])-1\n",
    "\n",
    "dir_ids = [dir_id_1, dir_id_2, dir_id_3]\n",
    "neu_ids = [neu_id_1]\n",
    "\n",
    "dir_edges = set([tuple(np.array(edge['node-number-list'])-1) for marker in dir_ids for edge in boundaryElements[marker]])\n",
    "neu_edges = set([tuple(np.array(edge['node-number-list'])-1) for marker in neu_ids for edge in boundaryElements[marker]])\n",
    "\n",
    "######################\n",
    "#                    #\n",
    "#     Hole edges     #\n",
    "#                    #\n",
    "######################\n",
    "# first hole\n",
    "hole_neu_edge_1 = np.array([edge['node-number-list'] for edge in boundaryElements[hole_neu_id_1]])-1\n",
    "hole_neu_edge_2 = np.array([edge['node-number-list'] for edge in boundaryElements[hole_neu_id_2]])-1\n",
    "hole_neu_edge_3 = np.array([edge['node-number-list'] for edge in boundaryElements[hole_neu_id_3]])-1\n",
    "hole_neu_edge_4 = np.array([edge['node-number-list'] for edge in boundaryElements[hole_neu_id_4]])-1\n",
    "\n",
    "hole_neu_edges = [hole_neu_edge_1, hole_neu_edge_2, hole_neu_edge_3, hole_neu_edge_4]\n",
    "\n",
    "ccw_hole_neu_edges = np.array([edge[::-1] for edges in hole_neu_edges[::-1] for edge in edges[::-1]])\n",
    "\n",
    "ccw_hole_neu_edge_1 = np.array([edge[::-1] for edge in hole_neu_edge_1[::-1]])\n",
    "ccw_hole_neu_edge_2 = np.array([edge[::-1] for edge in hole_neu_edge_2[::-1]])\n",
    "ccw_hole_neu_edge_3 = np.array([edge[::-1] for edge in hole_neu_edge_3[::-1]])\n",
    "ccw_hole_neu_edge_4 = np.array([edge[::-1] for edge in hole_neu_edge_4[::-1]])\n",
    "\n",
    "print(ccw_hole_neu_edges)\n",
    "print(ccw_hole_neu_edge_1)\n",
    "print(ccw_hole_neu_edge_2)\n",
    "print(ccw_hole_neu_edge_3)\n",
    "print(ccw_hole_neu_edge_4)\n",
    "\n",
    "#################################\n",
    "#                               #\n",
    "#  Crossing elements and edges  #\n",
    "#                               #\n",
    "#################################\n",
    "direction = np.array([0, -1])\n",
    "Bh = []\n",
    "Eh = {}\n",
    "\n",
    "Bh.append(tuple(hole_neu_edge_1[0]))\n",
    "\n",
    "# crossing flag \n",
    "connected = False\n",
    "\n",
    "######################################################\n",
    "# TODO: New structures of Eh, Bh for multiple holes. #\n",
    "######################################################\n",
    "while not connected:\n",
    "    cross_edge = Bh[-1]\n",
    "    # print(cross_edge)\n",
    "    # print(node[cross_edge[0]])\n",
    "    # print(node[cross_edge[1]])\n",
    "    if (cross_edge not in dir_edges) and (cross_edge not in neu_edges):\n",
    "        cross_edge = cross_edge[::-1]\n",
    "        ei = edge[cross_edge]\n",
    "        \n",
    "        i, j, k = element[ei]\n",
    "        order = [i,j,k]  \n",
    "        \n",
    "        local_edges = set([(i,j), (j,k), (k,i)])\n",
    "        local_edges.remove(cross_edge)\n",
    "        e1 = np.array(local_edges.pop())\n",
    "        e2 = np.array(local_edges.pop())\n",
    "        \n",
    "        e1n = node[e1[1]] - node[e1[0]]\n",
    "        e1n = e1n/np.linalg.norm(e1n)\n",
    "        \n",
    "        e2n = node[e2[1]] - node[e2[0]]\n",
    "        e2n = e2n/np.linalg.norm(e2n)\n",
    "        \n",
    "        d1 = np.abs(np.dot(direction, e1n))\n",
    "        d2 = np.abs(np.dot(direction, e2n))\n",
    "        if d1 < d2:\n",
    "            ae = np.zeros(3)\n",
    "            L1 = np.linalg.norm(node[cross_edge[0]] - node[cross_edge[1]])\n",
    "            L2 = np.linalg.norm(node[e1[0]] - node[e1[1]])\n",
    "            ae[order.index(cross_edge[0])]= 1/L1\n",
    "            ae[order.index(e1[0])]= -1/L2\n",
    "            Eh[ei] = ae\n",
    "            Bh.append(tuple(e1))\n",
    "        else:\n",
    "            ae = np.zeros(3)\n",
    "            L1 = np.linalg.norm(node[cross_edge[0]] - node[cross_edge[1]])\n",
    "            L2 = np.linalg.norm(node[e2[0]] - node[e2[1]])\n",
    "            ae[order.index(cross_edge[0])]= 1/L1\n",
    "            ae[order.index(e2[0])]= -1/L2\n",
    "            Eh[ei] = ae\n",
    "            Bh.append(tuple(e2))\n",
    "            \n",
    "        # print(e1,'\\t',e1n)\n",
    "        # print(e2,'\\t',e2n)\n",
    "        # print(ae)\n",
    "        # print(d1,d2)\n",
    "        # print(\"----------------------------------\")\n",
    "        \n",
    "            \n",
    "    else:\n",
    "        if cross_edge in neu_edges:\n",
    "            print('Neumann ! ')\n",
    "        else:\n",
    "            print('Dirichlet !')\n",
    "        connected = True\n",
    "    \n",
    "# print(Eh)\n",
    "        \n",
    "# np.all()\n",
    "\n",
    "######################\n",
    "#                    #\n",
    "#  Numbers of sets   #\n",
    "#                    #\n",
    "######################\n",
    "E = len(element)\n",
    "N = len(node)\n",
    "num_dof = N + nH\n",
    "\n",
    "# Display mesh\n",
    "cfv.draw_mesh(coords=coords, edof=edof, dofs_per_node=mesh.dofs_per_node,\n",
    "              el_type=mesh.el_type,filled=True,title=\"Example 2 - Mesh\")\n",
    "cfv.showAndWait()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "edf91826-efaf-416e-9306-d783aa7ea2c6",
   "metadata": {
    "tags": []
   },
   "source": [
    "## 2. Boundary condition"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "cb26e870-b29d-41b5-bb21-0c83e022121b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Externel boundary condition\n",
    "########## Dirichlet ##########\n",
    "def dir1(x, y):\n",
    "    return np.zeros_like(x)\n",
    "\n",
    "def dir2(x, y):\n",
    "    return np.sin(x)/np.sin(1)\n",
    "\n",
    "def dir3(x, y):\n",
    "    return np.zeros_like(x)\n",
    "########## Neumann ############\n",
    "def neu1(x, y):\n",
    "    return (1/np.tan(1))*np.sinh(y)/np.sinh(1)\n",
    "\n",
    "# Hole's boundary condition\n",
    "########## Dirichlet ##########\n",
    "# ...\n",
    "########## Neumann ############\n",
    "def hole_neu1(x, y):\n",
    "    return (np.sin(x)/np.sin(1))*(np.cosh(0.4)/np.sinh(1))\n",
    "\n",
    "def hole_neu2(x, y):\n",
    "    return -(np.cos(0.6)/np.sin(1))*(np.sinh(y)/np.sinh(1))\n",
    "\n",
    "def hole_neu3(x, y):\n",
    "    return -(np.sin(x)/np.sin(1))*(np.cosh(0.6)/np.sinh(1))\n",
    "\n",
    "def hole_neu4(x, y):\n",
    "    return (np.cos(0.4)/np.sin(1))*(np.sinh(y)/np.sinh(1))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d5fe1912-f8d3-4895-9e96-3156355eb71c",
   "metadata": {},
   "source": [
    "## 3. Assemble flexibility matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "fb18dbeb-89fc-4566-b9bc-7627f4bef33f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "size: 3361 x 3361\n",
      "rank: 3360\n"
     ]
    }
   ],
   "source": [
    "F = np.zeros((num_dof, num_dof))\n",
    "\n",
    "for i, e in enumerate(element):\n",
    "    # print(e)\n",
    "    # print(node.shape)\n",
    "    x1,y1 = node[e][0]\n",
    "    x2,y2 = node[e][1]\n",
    "    x3,y3 = node[e][2]\n",
    "    \n",
    "    A_e = 1/2*np.linalg.det(np.c_[node[e], np.ones((3,1))])\n",
    "    \n",
    "    Phi_e = (1/(2*A_e))*np.array([[x2-x3,x3-x1,x1-x2],\n",
    "                                  [y2-y3,y3-y1,y1-y2]])\n",
    "    \n",
    "    \n",
    "    if i in Eh:\n",
    "        L1 = np.sqrt((x2-x1)**2+(y2-y1)**2)\n",
    "        L2 = np.sqrt((x3-x2)**2+(y3-y2)**2)\n",
    "\n",
    "        left = np.linalg.inv([[(y1-y2)/L1,(x2-x1)/L1],\n",
    "                              [(y2-y3)/L2,(x3-x2)/L2]])\n",
    "\n",
    "        Phi_L = np.c_[left, np.zeros((2,1))]\n",
    "        \n",
    "        # print(Eh[i])\n",
    "        Phi_e = np.c_[Phi_e, Phi_L @ Eh[i]]\n",
    "        F_e = A_e * Phi_e.T @ Phi_e\n",
    "        \n",
    "        ### Todo: multiple holes ###\n",
    "        ey = np.append(e, -1)\n",
    "\n",
    "        F[np.ix_(ey,ey)] = F[np.ix_(ey,ey)] + F_e\n",
    "    else:\n",
    "        F_e = A_e * Phi_e.T @ Phi_e\n",
    "        \n",
    "        F[np.ix_(e,e)] = F[np.ix_(e,e)] + F_e\n",
    "\n",
    "print('size: {} x {}'.format(*F.shape))\n",
    "print('rank: {}'.format(np.linalg.matrix_rank(F)))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7a01a3f7-c7e0-40b2-aeef-6ac3d2787a52",
   "metadata": {},
   "source": [
    "## 4. Assemble kinematic vector (with source, Dirichlet problem)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "c855dbeb-fe00-4b49-8f23-8514e2fd3e53",
   "metadata": {},
   "outputs": [],
   "source": [
    "du = np.zeros(num_dof)\n",
    "\n",
    "##########################\n",
    "#   Dirichlet boundary   #\n",
    "##########################\n",
    "for edge in dir_edge_1:\n",
    "    i, j = edge\n",
    "    L = np.linalg.norm(node[i]-node[j])\n",
    "    up = gauss_line_int(dir1, *node[i], *node[j])\n",
    "    du[i] = du[i]-up/L\n",
    "    du[j] = du[j]+up/L\n",
    "\n",
    "for edge in dir_edge_2:\n",
    "    i, j = edge\n",
    "    L = np.linalg.norm(node[i]-node[j])\n",
    "    up = gauss_line_int(dir2, *node[i], *node[j])\n",
    "    du[i] = du[i]-up/L\n",
    "    du[j] = du[j]+up/L\n",
    "\n",
    "for edge in dir_edge_3:\n",
    "    i, j = edge\n",
    "    L = np.linalg.norm(node[i]-node[j])\n",
    "    up = gauss_line_int(dir3, *node[i], *node[j])\n",
    "    du[i] = du[i]-up/L\n",
    "    du[j] = du[j]+up/L\n",
    "\n",
    "##########################################\n",
    "#    Neumann boundary (without source)   #\n",
    "##########################################\n",
    "s_N = np.zeros(num_dof)\n",
    "A = np.eye(num_dof)\n",
    "\n",
    "num_red = len(neu_edge_1) + len(ccw_hole_neu_edges)\n",
    "num_eff = num_dof - num_red\n",
    "\n",
    "# 1. Externel edge\n",
    "# effective variable (s1) is at first node of first Neumann edge \n",
    "# (edge's nodes are ordered counterclockwise within the local element).\n",
    "ext_eff_dofs = [neu_edge_1[0][0]]\n",
    "ext_red_dofs = [edge[1] for edge in neu_edge_1] \n",
    "\n",
    "# 2. Hole edge\n",
    "# last edge should be excluded for closed Neumann boundary. \n",
    "hole_eff_dofs = [hole_neu_edge_1[0][0]]\n",
    "hole_red_dofs = [edge[1] for edge in ccw_hole_neu_edges[:-1]]\n",
    "\n",
    "red_dofs = ext_red_dofs + hole_red_dofs\n",
    "\n",
    "A[np.ix_(ext_red_dofs, ext_eff_dofs)] = 1\n",
    "A[np.ix_(hole_red_dofs, hole_eff_dofs)] = 1\n",
    "\n",
    "A = np.delete(A, red_dofs, axis=1)\n",
    "\n",
    "# Neumann terms\n",
    "####### TODO : crossing edge #########\n",
    "# externel (counter-clockwise)\n",
    "fl_bar = 0\n",
    "for edge in neu_edge_1:\n",
    "    i, j = edge\n",
    "    fl_bar = fl_bar + gauss_line_int(neu1, *node[i], *node[j])\n",
    "    s_N[j] = fl_bar\n",
    "\n",
    "# hole (clockwise)\n",
    "# last edge should be excluded for closed Neumann boundary. \n",
    "fl_bar = 0\n",
    "for edge in ccw_hole_neu_edge_4:\n",
    "    i, j = edge\n",
    "    fl_bar = fl_bar + gauss_line_int(hole_neu4, *node[i], *node[j])\n",
    "    s_N[j] = fl_bar\n",
    "    \n",
    "for edge in ccw_hole_neu_edge_3:\n",
    "    i, j = edge\n",
    "    fl_bar = fl_bar + gauss_line_int(hole_neu3, *node[i], *node[j])\n",
    "    s_N[j] = fl_bar\n",
    "    \n",
    "for edge in ccw_hole_neu_edge_2:\n",
    "    i, j = edge\n",
    "    fl_bar = fl_bar + gauss_line_int(hole_neu2, *node[i], *node[j])\n",
    "    s_N[j] = fl_bar\n",
    "    \n",
    "for edge in ccw_hole_neu_edge_1[:-1]:\n",
    "    i, j = edge\n",
    "    fl_bar = fl_bar + gauss_line_int(hole_neu1, *node[i], *node[j])\n",
    "    s_N[j] = fl_bar\n",
    "    \n",
    "# Source \n",
    "# Particular solution q_0\n",
    "# def qx_0(x, y):\n",
    "#     return x\n",
    "\n",
    "# def qy_0(x, y):\n",
    "#     return np.zeros_like(x)\n",
    "\n",
    "# for i, e in enumerate(element):\n",
    "#     x1,y1 = node[e][0]\n",
    "#     x2,y2 = node[e][1]\n",
    "#     x3,y3 = node[e][2]\n",
    "    \n",
    "#     A_e = 1/2*np.linalg.det(np.c_[node[e], np.ones((3,1))])\n",
    "    \n",
    "#     Phi_e = (1/(2*A_e))*np.array([[x2-x3,x3-x1,x1-x2],\n",
    "#                                   [y2-y3,y3-y1,y1-y2]])\n",
    "    \n",
    "#     int_qx = gauss_tri_int(qx_0, x1, y1, x2, y2, x3, y3, degree=2)\n",
    "#     int_qy = gauss_tri_int(qy_0, x1, y1, x2, y2, x3, y3, degree=2)\n",
    "    \n",
    "#     par = Phi_e.T @ np.array([int_qx, int_qy])\n",
    "    \n",
    "#     du[e] = du[e] - par\n",
    "\n",
    "# Neumann boundary with source.\n",
    "# ..."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5619f835-d960-4905-a43c-5b77ff0dcc40",
   "metadata": {},
   "source": [
    "## 5. Constraint"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "0dba5f60-fa70-4d41-ac50-72d386d3aa69",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3.469446951953614e-18\n"
     ]
    }
   ],
   "source": [
    "FU = F[:]\n",
    "dU = du[:]\n",
    "\n",
    "F = A.T @ FU @ A\n",
    "d = A.T @ (dU - FU@s_N)\n",
    "\n",
    "r1 = 0\n",
    "\n",
    "for edge in hole_neu_edge_1:\n",
    "    i, j = edge\n",
    "    r1 = r1 + gauss_line_int(hole_neu1, *node[i], *node[j])\n",
    "\n",
    "for edge in hole_neu_edge_2:\n",
    "    i, j = edge\n",
    "    r1 = r1 + gauss_line_int(hole_neu2, *node[i], *node[j])\n",
    "\n",
    "for edge in hole_neu_edge_3:\n",
    "    i, j = edge\n",
    "    r1 = r1 + gauss_line_int(hole_neu3, *node[i], *node[j])\n",
    "\n",
    "for edge in hole_neu_edge_4:\n",
    "    i, j = edge\n",
    "    r1 = r1 + gauss_line_int(hole_neu4, *node[i], *node[j])\n",
    "    \n",
    "print(r1)\n",
    "\n",
    "F[-1,:] = 0\n",
    "F[-1,-1] = 1\n",
    "d[-1] = r1\n",
    "\n",
    "F[0,:] = 0\n",
    "F[0,0] = 1\n",
    "d[0] = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "2c04f820-d63c-4dbd-9f7d-00b07d8bec82",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x2be2cacd670>"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkEAAAI/CAYAAABwLA0cAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAeUklEQVR4nO3dbWyd53nY8euyaeOQVIFENENZjDN2SrAmHVZ7EIwKNYasXYu0X5ICRdEACzJYgPWhAVKsHxb0S9NhA7qhL/syZHZBoV7RNg2adA2GYGsQBIgKSGmtyM2b2TkOkjiSRcnUgsYiVMf2vQ866VRVsijeD3kOdf1+gCDy4bnvc+XJsfTHc16UrbUAAKjmrkkPAAAwCSIIAChJBAEAJYkgAKAkEQQAlCSCAICSZnbzzu699942Go22vf573/vegNMAABVcuXLlxdba4vXHdzWCRqNRPPzww9tef/bs2QGnAQAqWFtb++aNjns6DAAoSQQBACV1RVBmvisz/zozv5aZHxpqKACAnbbtCMrMuyPiv0bET0fEOyLivZn5jqEGAwDYST1Xgh6OiK+11r7eWns5Ij4aEe8eZiwAgJ3VE0HLEfH8Nd9/e3wMAGDq7fgLozPzscx8KjOf8jk/AMC06ImgsxHxwDXfv3l87O9prT3RWjvcWjt8zz33dNwdAMBweiLoLyPibZn5g5l5b0T8QkR8cpixAAB21rY/Mbq19kpmfiAi/ndE3B0Rx1trXxlsMgCAHdT1z2a01j4VEZ8aaBYAgF3jE6MBgJJEEABQkggCAErqek3Q7fre974XZ8/+g3fRb9nv/d7vdc/wvve9r3sPIq5cuTLpEWI0Gk16BAb0/PPP3/pGr+OBBx649Y3YFS+++GL3Hvfdd1/X+m9961vdM7zlLW/p3qPXs88+27U+M7vWv/Wtb+1aP+1cCQIAShJBAEBJIggAKEkEAQAliSAAoCQRBACUJIIAgJJEEABQkggCAEoSQQBASSIIAChJBAEAJYkgAKAkEQQAlCSCAICSRBAAUFK21nbtzmZnZ9vKysqu3d+NvOlNb+paf+HChYEm4U5w5cqVSY8Qo9Fo0iMA3NTly5e71s/Pz3fPsLa2drq1dvj6464EAQAliSAAoCQRBACUJIIAgJJEEABQkggCAEoSQQBASSIIAChJBAEAJYkgAKAkEQQAlCSCAICSRBAAUJIIAgBKEkEAQEkzkx5gt124cKFr/TPPPNM9w9vf/vbuPZgOo9Fo0iPcMTY3Nyc9wlR44YUXJj3CVLj//vsnPcJUmJubm/QI3ebn57vWv+ENbxhmkBtwJQgAKEkEAQAliSAAoCQRBACUJIIAgJJEEABQkggCAEoSQQBASSIIAChJBAEAJYkgAKAkEQQAlCSCAICSRBAAUJIIAgBKEkEAQEkzu3lnr732Wly5cmXb60ej0YDTbM/b3/727j0eeuihrvUnT57snmEaziVca25ubtIjTIVDhw5NegQGsrm5OekRup09e7Z7j+Xl5a71Pd1wK64EAQAliSAAoCQRBACUJIIAgJJEEABQkggCAEoSQQBASSIIAChJBAEAJYkgAKAkEQQAlCSCAICSRBAAUJIIAgBKEkEAQEkzu3lnd911V4xGo928y6l08uTJrvVHjhzpnuHEiRNd6/ft29c9AzCdNjc3u9bPzc0NNMnedvny5e49Jn0ul5eXJ3r/ERHPPvvsju3tShAAUJIIAgBKEkEAQEkiCAAoqeuF0Zn5jYj4bkS8GhGvtNYODzEUAMBOG+LdYf+ytfbiAPsAAOwaT4cBACX1RlCLiD/LzNOZ+dgQAwEA7Ibep8Meaa2dzcw3RcSnM3Ottfa5a28wjqPHIiJmZnb1sxkBAG6q60pQa+3s+PcLEfEnEfHwDW7zRGvtcGvtsAgCAKbFtiMoM+cz8we+/3VE/FREfHmowQAAdlLPpZmliPiTzPz+Pn/QWvtfg0wFALDDth1BrbWvR8SPDDgLAMCu8RZ5AKAkEQQAlCSCAICSsrW2a3c2OzvbVlZWdu3+7lQvvfRS9x6PPPJI1/pTp051zzAajbr3gO/b2Njo3mNhYWGASfa+y5cvd62fn58faBIYxtra2ukb/fumrgQBACWJIACgJBEEAJQkggCAkkQQAFCSCAIAShJBAEBJIggAKEkEAQAliSAAoCQRBACUJIIAgJJEEABQkggCAEoSQQBASTOTHuB2XLlypXuP0Wg0wCSTtW/fvu49Tp061bX+yJEj3TOcOXOmew/4vjvhv+1pMT8/P+kRGDt37lzX+oMHDw40yZ3JlSAAoCQRBACUJIIAgJJEEABQkggCAEoSQQBASSIIAChJBAEAJYkgAKAkEQQAlCSCAICSRBAAUJIIAgBKEkEAQEkiCAAoSQQBACXNTHoAJmM0GnWtP3PmTPcMDz30UNf6kydPds/Qex6YHvPz85MegQFtbm52rZ+bmxtoksk6ePBg1/qzZ892rV9eXu5aP+1cCQIAShJBAEBJIggAKEkEAQAliSAAoCQRBACUJIIAgJJEEABQkggCAEoSQQBASSIIAChJBAEAJYkgAKAkEQQAlCSCAICSZiY9wO0YjUaTHoEBnTx5smv9kSNHumc4ceJE1/p9+/Z1zwBD2tzc7N5jbm5ugEn6XL58uWv9NPxvmAbLy8uTHmGquRIEAJQkggCAkkQQAFCSCAIAShJBAEBJIggAKEkEAQAliSAAoCQRBACUJIIAgJJEEABQkggCAEoSQQBASSIIAChJBAEAJYkgAKCkmUkPQF2j0ahr/YkTJ7pneOSRR7rWnzp1qnuG3vPA9NjY2OjeY2FhoWt9a617hmmwuLg46REowJUgAKAkEQQAlCSCAICSbhlBmXk8My9k5pevObY/Mz+dmc+Of3/jzo4JADCsrVwJ+t2IeNd1xz4UEZ9prb0tIj4z/h4AYM+4ZQS11j4XEZeuO/zuiHhy/PWTEfGeYccCANhZ231N0FJr7YXx1+cjYmmgeQAAdkX35wS11lpm3vSDKTLzsYh4LCJiZsbHEgEA02G7V4LWM/P+iIjx7xdudsPW2hOttcOttcMiCACYFtuNoE9GxPvHX78/Iv50mHEAAHbHVt4i/4cRcTIi/klmfjszj0bEr0fET2bmsxHxr8bfAwDsGbd8fqq19t6b/OgnBp4FAGDX+MRoAKAkEQQAlCSCAICSyr1n/cqVK13rR6PRQJPQa9++fd17nDp1qmv9kSNHumc4c+ZM9x5Mh2n482F+fn7SI8Dfc+7cua71Bw8eHGiSf8iVIACgJBEEAJQkggCAkkQQAFCSCAIAShJBAEBJIggAKEkEAQAliSAAoCQRBACUJIIAgJJEEABQkggCAEoSQQBASSIIAChJBAEAJc1MegCYpNFo1LX+zJkz3TM89NBDXetPnjzZPUPveeCq+fn5SY/AgDY3N7vWz83NDTTJ3tZam/QIN+VKEABQkggCAEoSQQBASSIIAChJBAEAJYkgAKAkEQQAlCSCAICSRBAAUJIIAgBKEkEAQEkiCAAoSQQBACWJIACgJBEEAJQ0M+kBdttoNJr0CPD3nDx5smv9kSNHumc4ceJE1/p9+/Z1z8BVm5ubXevn5uYGmoTLly93rff/xVXLy8uTHuGmXAkCAEoSQQBASSIIAChJBAEAJYkgAKAkEQQAlCSCAICSRBAAUJIIAgBKEkEAQEkiCAAoSQQBACWJIACgJBEEAJQkggCAkkQQAFDSzKQHgOpGo1HX+hMnTnTP8Mgjj3StP3XqVPcMvedhCBsbG13rFxYWumdorXXvwTAWFxcnPQI7zJUgAKAkEQQAlCSCAICSRBAAUJIIAgBKEkEAQEkiCAAoSQQBACWJIACgJBEEAJQkggCAkkQQAFCSCAIAShJBAEBJIggAKGlmN+/stddeiytXrmx7/Wg0GnAauDPs27eve49Tp051rT9y5Ej3DGfOnOneo9c0/BkzPz8/6RHg75w7d657j4MHDw4wyc5wJQgAKEkEAQAliSAAoCQRBACUdMsIyszjmXkhM798zbEPZ+bZzHx6/OtndnZMAIBhbeVK0O9GxLtucPy3W2sPjn99atixAAB21i0jqLX2uYi4tAuzAADsmp7XBH0gM784frrsjYNNBACwC7YbQR+JiEMR8WBEvBARv3mzG2bmY5n5VGY+9eqrr27z7gAAhrWtCGqtrbfWXm2tvRYRvxMRD7/ObZ9orR1urR2+++67tzsnAMCgthVBmXn/Nd/+bER8+Wa3BQCYRrf8t8My8w8j4p0RcV9mfjsifjUi3pmZD0ZEi4hvRMSxnRsRAGB4t4yg1tp7b3B4dQdmAQDYNT4xGgAoSQQBACWJIACgpFu+JmhId911V4xGo22vf/7557tneOCBB7rWb25uds8wNzfXvQcMqee/y4iIM2fOdM9w9OjRrvWrq/0vVZyfn+/eg4iLFy9279H7kSr79+/vnuFOsL6+3rX+4MGDA00ynVwJAgBKEkEAQEkiCAAoSQQBACWJIACgJBEEAJQkggCAkkQQAFCSCAIAShJBAEBJIggAKEkEAQAliSAAoCQRBACUJIIAgJJmJj3A7XjggQcmPcIdY2Njo2v9aDTqnmF+fr57D+4cq6urXeuPHj068RnuFBcvXuxav7i4ONAk23fp0qXuPfbv3z/AJH3W19e71i8tLQ00yZ3JlSAAoCQRBACUJIIAgJJEEABQkggCAEoSQQBASSIIAChJBAEAJYkgAKAkEQQAlCSCAICSRBAAUJIIAgBKEkEAQEkiCAAoSQQBACVla23X7mx2dratrKzs2v3thM3Nze495ubmBpgEhrOxsdG1fjQadc8wPz/fvUevo0ePdq1fXV3tnuHixYtd6xcXF7tnYBiXLl3q3mP//v0DTLJ96+vr3XssLS0NMEmftbW10621w9cfdyUIAChJBAEAJYkgAKAkEQQAlCSCAICSRBAAUJIIAgBKEkEAQEkiCAAoSQQBACWJIACgJBEEAJQkggCAkkQQAFCSCAIASsrW2q7d2ezsbFtZWdm1+9sJm5ub3XvMzc0NMAlD2NjY6Fo/Go26Z5ifn+/eg+lw9OjR7j1WV1e71l+8eLF7hsXFxe497gSXLl3qWr9///6BJtnb1tfXu9YvLS11z7C2tna6tXb4+uOuBAEAJYkgAKAkEQQAlCSCAICSRBAAUJIIAgBKEkEAQEkiCAAoSQQBACWJIACgJBEEAJQkggCAkkQQAFCSCAIAShJBAEBJIggAKGlm0gPsNS+88EL3HocOHRpgEoawsLAw6RG6bW5udu/RWutaPz8/3z3DNLh48WLX+tXV1e4ZHn300a71x48f756Bq15++eVJj3BHWFpamvQIN+VKEABQkggCAEoSQQBASbeMoMx8IDM/m5lfzcyvZOYHx8f3Z+anM/PZ8e9v3PlxAQCGsZUrQa9ExC+31t4RET8aEb+Yme+IiA9FxGdaa2+LiM+MvwcA2BNuGUGttRdaa18Yf/3diHgmIpYj4t0R8eT4Zk9GxHt2aEYAgMHd1muCMnMlIh6KiM9HxFJr7fvvFz8fEdP7HjgAgOts+XOCMnNfRHw8In6ptfY3mfl3P2uttcy84QeNZOZjEfFYRMTMjI8lAgCmw5auBGXmPXE1gH6/tfaJ8eH1zLx//PP7I+LCjda21p5orR1urR0WQQDAtNjKu8MyIlYj4pnW2m9d86NPRsT7x1+/PyL+dPjxAAB2xlYuzfxYRLwvIr6UmU+Pj/1KRPx6RHwsM49GxDcj4ud3ZEIAgB1wywhqrf15RORNfvwTw44DALA7fGI0AFCSCAIAShJBAEBJ3rPOnrW5udm9R2s3/HirLZufn++eodfc3NykR7hjLC4udq2/ePFi9wzHjx/vWv/oo49OfIYhnD9/vmv9gQMHumcYYg+mmytBAEBJIggAKEkEAQAliSAAoCQRBACUJIIAgJJEEABQkggCAEoSQQBASSIIAChJBAEAJYkgAKAkEQQAlCSCAICSRBAAUJIIAgBKytbart3Z7OxsW1lZ2bX72wnPPfdc9x6HDh3qWr+5udk9Q+//7/Pz890zwJAuXrzYvcfi4uIAk+x9jz76aNf648ePDzTJ3nb+/PnuPQ4cODDAJHvb2traENucbq0dvv6gK0EAQEkiCAAoSQQBACWJIACgJBEEAJQkggCAkkQQAFCSCAIAShJBAEBJIggAKEkEAQAliSAAoCQRBACUJIIAgJJEEABQUrbWdu3OZmdn28rKyrbXv/jii90z3HfffV3rn3vuue4ZDh061L0H3GkuXrzYtX5xcXGgSej16KOPdu9x/PjxrvXnz5/vnuHAgQPdezAd1tbWTrfWDl9/3JUgAKAkEQQAlCSCAICSRBAAUJIIAgBKEkEAQEkiCAAoSQQBACWJIACgJBEEAJQkggCAkkQQAFCSCAIAShJBAEBJIggAKEkEAQAlzUx6gNtx3333TXqEuP/++yc9wh1jc3Oza/3ly5e7Z1hcXOzeg2Hcfffdkx7hjnD+/PnuPQ4cONC1/vjx490zHD16tGv96upq9wx3go2Nja71CwsLA00ynVwJAgBKEkEAQEkiCAAoSQQBACWJIACgJBEEAJQkggCAkkQQAFCSCAIAShJBAEBJIggAKEkEAQAliSAAoCQRBACUJIIAgJJmJj3A7fjWt77Vvcdb3vKWASZhCHNzcxNdz3TZv39/1/pLly5NfIZpcODAgUmPEOfPn+/eY3V1tWv90aNHJz7DEDY2NrrWLywsDDTJncmVIACgJBEEAJQkggCAkm4ZQZn5QGZ+NjO/mplfycwPjo9/ODPPZubT418/s/PjAgAMYysvjH4lIn65tfaFzPyBiDidmZ8e/+y3W2u/sXPjAQDsjFtGUGvthYh4Yfz1dzPzmYhY3unBAAB20m29JigzVyLioYj4/PjQBzLzi5l5PDPfOPRwAAA7ZcsRlJn7IuLjEfFLrbW/iYiPRMShiHgwrl4p+s2brHssM5/KzKdeeeWV/okBAAawpQjKzHviagD9fmvtExERrbX11tqrrbXXIuJ3IuLhG61trT3RWjvcWjs8M7OnPpsRALiDbeXdYRkRqxHxTGvtt645fv81N/vZiPjy8OMBAOyMrVya+bGIeF9EfCkznx4f+5WIeG9mPhgRLSK+ERHHdmA+AIAdsZV3h/15ROQNfvSp4ccBANgdPjEaAChJBAEAJYkgAKAkEQQAlJSttV27s9nZ2baysrJr97cTNjc3u/eYm5sbYBLuFL2PqcuXL3fPsLi42L0HEZcuXere4+WXX+5af+DAge4ZuOro0aNd61dXVweaZG/b2NjoWr+wsNA9w9ra2unW2uHrj7sSBACUJIIAgJJEEABQkggCAEoSQQBASSIIAChJBAEAJYkgAKAkEQQAlCSCAICSRBAAUJIIAgBKEkEAQEkiCAAoSQQBACXNTHqAijY3N7vWX758uXuGxcXF7j0Yxtzc3ETX8/9dunSpa/3+/fsHmmSyzp8/37X+wIEDA00yWaurq13rjx49OvEZem1sbHTvsbCw0LX+h3/4h7tnWFtbu+FxV4IAgJJEEABQkggCAEoSQQBASSIIAChJBAEAJYkgAKAkEQQAlCSCAICSRBAAUJIIAgBKEkEAQEkiCAAoSQQBACWJIACgJBEEAJSUrbVdu7PZ2dm2srKy7fXPPvts9wxve9vbuvcg4ty5c917HDx4cIBJ4Kr19fXuPZaWlgaYhGmwsbHRvcfCwsIAk/Q5duxY1/rHH398oEn2trW1tdOttcPXH3clCAAoSQQBACWJIACgJBEEAJQkggCAkkQQAFCSCAIAShJBAEBJIggAKEkEAQAliSAAoCQRBACUJIIAgJJEEABQkggCAEqamfQAtyMzJz0CYwcPHuze4+zZs13rl5eXu2dgeqyvr3etX1paGmgSpsHGxkbX+oWFhYEmmazHH3+8a/2xY8cmev/TzpUgAKAkEQQAlCSCAICSRBAAUJIIAgBKEkEAQEkiCAAoSQQBACWJIACgJBEEAJQkggCAkkQQAFCSCAIAShJBAEBJIggAKEkEAQAlZWtt1+7s3nvvbQcOHNj2+vn5+QGnmZyzZ892rV9eXh5oEhjGuXPnuvc4ePDgAJP0WV9f71q/tLQ00CTcCTY2Nrr3WFhYGGCS7Tt27Fj3Ho8//vgAk/RZW1s73Vo7fP1xV4IAgJJEEABQkggCAEq6ZQRl5igz/yIz/yozv5KZvzY+/oOZ+fnM/Fpm/lFm3rvz4wIADGMrV4L+NiJ+vLX2IxHxYES8KzN/NCL+U0T8dmvtrRHxfyPi6I5NCQAwsFtGULvqpfG394x/tYj48Yj44/HxJyPiPTsxIADATtjSa4Iy8+7MfDoiLkTEpyPiuYj4TmvtlfFNvh0R3rcNAOwZW4qg1tqrrbUHI+LNEfFwRPzQVu8gMx/LzKcy86nXXntte1MCAAzstt4d1lr7TkR8NiKORMQbMnNm/KM3R8QNPwGwtfZEa+1wa+3wXXd5MxoAMB228u6wxcx8w/jr2Yj4yYh4Jq7G0M+Nb/b+iPjTHZoRAGBwM7e+SdwfEU9m5t1xNZo+1lr7n5n51Yj4aGb+h4g4ExGrOzgnAMCgbhlBrbUvRsRDNzj+9bj6+iAAgD3Hi3QAgJJEEABQkggCAErK1tqu3dns7GxbWVnZtfu7kcuXL3etn5+fH2gS7gTnzp3r3qP3v8HlZZ9TOi3W19e791haWhpgEoawsbHRtX5hYWGgSfa2Y8eOda1//PHHu2dYW1s73Vo7fP1xV4IAgJJEEABQkggCAEoSQQBASSIIAChJBAEAJYkgAKAkEQQAlCSCAICSRBAAUJIIAgBKEkEAQEkiCAAoSQQBACWJIACgpGyt7dqdzc7OtpWVlV27PwBgb1teXu7e4zOf+czp1trh64+7EgQAlCSCAICSRBAAUJIIAgBKEkEAQEkiCAAoSQQBACWJIACgJBEEAJQkggCAkkQQAFCSCAIAShJBAEBJIggAKEkEAQAliSAAoKRsre3enWVejIhvvs5N7ouIF3dpnDudczkM53EYzuNwnMthOI/D2Cvn8R+11havP7irEXQrmflUa+3wpOe4EziXw3Aeh+E8Dse5HIbzOIy9fh49HQYAlCSCAICSpi2Cnpj0AHcQ53IYzuMwnMfhOJfDcB6HsafP41S9JggAYLdM25UgAIBdMTURlJnvysy/zsyvZeaHJj3PXpWZ38jML2Xm05n51KTn2Usy83hmXsjML19zbH9mfjoznx3//sZJzrgX3OQ8fjgzz44fl09n5s9Mcsa9IDMfyMzPZuZXM/MrmfnB8XGPydvwOufRY/I2ZeYoM/8iM/9qfC5/bXz8BzPz8+O/v/8oM++d9KxbNRVPh2Xm3RHxfyLiJyPi2xHxlxHx3tbaVyc62B6Umd+IiMOttb3wuQ1TJTP/RUS8FBH/vbX2T8fH/nNEXGqt/fo4zt/YWvt3k5xz2t3kPH44Il5qrf3GJGfbSzLz/oi4v7X2hcz8gYg4HRHviYh/Ex6TW/Y65/Hnw2PytmRmRsR8a+2lzLwnIv48Ij4YEf82Ij7RWvtoZv63iPir1tpHJjnrVk3LlaCHI+JrrbWvt9ZejoiPRsS7JzwTxbTWPhcRl647/O6IeHL89ZNx9Q9PXsdNziO3qbX2QmvtC+OvvxsRz0TEcnhM3pbXOY/cpnbVS+Nv7xn/ahHx4xHxx+Pje+oxOS0RtBwRz1/z/bfDg3S7WkT8WWaezszHJj3MHWCptfbC+OvzEbE0yWH2uA9k5hfHT5d5Cuc2ZOZKRDwUEZ8Pj8ltu+48RnhM3rbMvDszn46ICxHx6Yh4LiK+01p7ZXyTPfX397REEMN5pLX2zyPipyPiF8dPTTCAdvW548k/f7w3fSQiDkXEgxHxQkT85kSn2UMyc19EfDwifqm19jfX/sxjcutucB49JrehtfZqa+3BiHhzXH0W54cmO1GfaYmgsxHxwDXfv3l8jNvUWjs7/v1CRPxJXH2Qsn3r49cUfP+1BRcmPM+e1FpbH//h+VpE/E54XG7J+HUXH4+I32+tfWJ82GPyNt3oPHpM9mmtfSciPhsRRyLiDZk5M/7Rnvr7e1oi6C8j4m3jV5jfGxG/EBGfnPBMe05mzo9f+BeZOR8RPxURX379VdzCJyPi/eOv3x8RfzrBWfas7/+lPfaz4XF5S+MXoa5GxDOttd+65kcek7fhZufRY/L2ZeZiZr5h/PVsXH0z0zNxNYZ+bnyzPfWYnIp3h0VEjN+e+F8i4u6ION5a+4+TnWjvycx/HFev/kREzETEHziPW5eZfxgR74yr/yryekT8akT8j4j4WES8JSK+GRE/31rzot/XcZPz+M64+rRDi4hvRMSxa17Xwg1k5iMRcSIivhQRr40P/0pcfT2Lx+QWvc55fG94TN6WzPxncfWFz3fH1YsoH2ut/fvx3z0fjYj9EXEmIv51a+1vJzfp1k1NBAEA7KZpeToMAGBXiSAAoCQRBACUJIIAgJJEEABQkggCAEoSQQBASSIIACjp/wEvTwBtjqXt4QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x720 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "plt.figure(figsize=(10,10))\n",
    "plt.imshow(F,cmap='gray')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4c4ffc26-a9f0-408c-80eb-723329c3b8ee",
   "metadata": {},
   "source": [
    "## 6. Solve and show"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "975aa11a-a4b7-49fb-869b-39843c1242e6",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3.469446951953614e-18\n",
      "3.469446951953614e-18\n",
      "1.4574548369358473e-12\n"
     ]
    }
   ],
   "source": [
    "s_eff = np.linalg.inv(F) @ d\n",
    "s = A @ s_eff + s_N\n",
    "print(s_eff[-1])\n",
    "print(d[-1])\n",
    "print( np.linalg.norm(F @ s_eff - d))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c8e4ab57-f7ab-4b05-a27f-61405d7d3079",
   "metadata": {},
   "source": [
    "## 7. Construct the flux field"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "d4a4a3c5-081a-43ff-9f2a-1603428996bf",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-0.00990582 -0.40433158]\n",
      "[-0.00990582 -0.40433158]\n",
      "[-0.35485114 -0.43218963]\n",
      "[-0.35485114 -0.43218963]\n",
      "[-0.05818146 -0.41133421]\n",
      "[-0.05818146 -0.41133421]\n",
      "[-0.0683943  -0.40546168]\n",
      "[-0.0683943  -0.40546168]\n",
      "[-0.32929458 -0.42642685]\n",
      "[-0.32929458 -0.42642685]\n",
      "[-0.33877401 -0.43197881]\n",
      "[-0.33877401 -0.43197881]\n",
      "[-0.36311161 -0.42788121]\n",
      "[-0.36311161 -0.42788121]\n",
      "[-0.32256208 -0.43026051]\n",
      "[-0.32256208 -0.43026051]\n",
      "[-0.26063426 -0.41741306]\n",
      "[-0.26063426 -0.41741306]\n",
      "[-0.27089501 -0.4234165 ]\n",
      "[-0.27089501 -0.4234165 ]\n",
      "[-0.04981387 -0.40606869]\n",
      "[-0.04981387 -0.40606869]\n",
      "[-0.31202442 -0.42410931]\n",
      "[-0.31202442 -0.42410931]\n",
      "[-0.27779171 -0.41948634]\n",
      "[-0.27779171 -0.41948634]\n",
      "[-0.24342794 -0.41554383]\n",
      "[-0.24342794 -0.41554383]\n",
      "[-0.20873802 -0.4119911 ]\n",
      "[-0.20873802 -0.4119911 ]\n",
      "[-0.25348755 -0.42147977]\n",
      "[-0.25348755 -0.42147977]\n",
      "[-0.0401875  -0.41175946]\n",
      "[-0.0401875  -0.41175946]\n",
      "[-0.13072967 -0.41212663]\n",
      "[-0.13072967 -0.41212663]\n",
      "[-0.21869324 -0.41796208]\n",
      "[-0.21869324 -0.41796208]\n",
      "[-0.20139975 -0.41636814]\n",
      "[-0.20139975 -0.41636814]\n",
      "[-0.14832414 -0.41284355]\n",
      "[-0.14832414 -0.41284355]\n",
      "[-0.30536229 -0.42792945]\n",
      "[-0.30536229 -0.42792945]\n",
      "[-0.28835201 -0.42563299]\n",
      "[-0.28835201 -0.42563299]\n",
      "[-0.23607585 -0.41972646]\n",
      "[-0.23607585 -0.41972646]\n",
      "[-0.10355267 -0.40567493]\n",
      "[-0.10355267 -0.40567493]\n",
      "[-0.07659503 -0.41072364]\n",
      "[-0.07659503 -0.41072364]\n",
      "[-0.12115727 -0.40652076]\n",
      "[-0.12115727 -0.40652076]\n",
      "[-0.13912253 -0.40739091]\n",
      "[-0.13912253 -0.40739091]\n",
      "[-0.0306515  -0.40607003]\n",
      "[-0.0306515  -0.40607003]\n",
      "[-0.22602776 -0.41376077]\n",
      "[-0.22602776 -0.41376077]\n",
      "[-0.29492795 -0.42181378]\n",
      "[-0.29492795 -0.42181378]\n",
      "[-0.11306458 -0.4112126 ]\n",
      "[-0.11306458 -0.4112126 ]\n",
      "[-0.17476188 -0.40874331]\n",
      "[-0.17476188 -0.40874331]\n",
      "[-0.18396846 -0.41476739]\n",
      "[-0.18396846 -0.41476739]\n",
      "[-0.15688032 -0.40800267]\n",
      "[-0.15688032 -0.40800267]\n",
      "[-0.08633427 -0.40535134]\n",
      "[-0.08633427 -0.40535134]\n",
      "[-0.19143242 -0.41008731]\n",
      "[-0.19143242 -0.41008731]\n",
      "[-0.09497866 -0.41066898]\n",
      "[-0.09497866 -0.41066898]\n",
      "[-0.0104066  -0.40369098]\n",
      "[-0.0104066  -0.40369098]\n",
      "[-0.1660419  -0.41375971]\n",
      "[-0.1660419  -0.41375971]\n",
      "[-0.37348518 -0.43297374]\n",
      "[-0.37348518 -0.43297374]\n",
      "[-0.37234496 -0.43408687]\n",
      "[-0.37234496 -0.43408687]\n",
      "[-0.02250151 -0.41144927]\n",
      "[-0.02250151 -0.41144927]\n",
      "[-0.34742472 -0.42765901]\n",
      "[-0.34742472 -0.42765901]\n",
      "<visvis.wibjects.colorWibjects.Colorbar object at 0x000002BE2F0896A0>\n",
      "ERROR calling '_OnAxesPositionChange':\n",
      "  File \"D:\\20220326(EFEM-Poisson)\\EFEM\\EFEM\\lib\\site-packages\\visvis\\wibjects\\colorWibjects.py\", line 622, in _OnAxesPositionChange\n",
      "    self.position.x = axes.position.width+5\n",
      "  File \"D:\\20220326(EFEM-Poisson)\\EFEM\\EFEM\\lib\\site-packages\\visvis\\core\\misc.py\", line 217, in fsetWithDraw\n",
      "    fset(self, *args)\n",
      "  File \"D:\\20220326(EFEM-Poisson)\\EFEM\\EFEM\\lib\\site-packages\\visvis\\core\\base.py\", line 1126, in fset\n",
      "    self._Update()\n",
      "  File \"D:\\20220326(EFEM-Poisson)\\EFEM\\EFEM\\lib\\site-packages\\visvis\\core\\base.py\", line 928, in _Update\n",
      "    self._CalculateInPixels()\n",
      "  File \"D:\\20220326(EFEM-Poisson)\\EFEM\\EFEM\\lib\\site-packages\\visvis\\core\\base.py\", line 997, in _CalculateInPixels\n",
      "    raise Exception(\"Can only calculate the position in pixels\"+\n",
      "Exception: Can only calculate the position in pixels if the owner has a parent!\n",
      "\n"
     ]
    }
   ],
   "source": [
    "q_h = np.zeros((E,2)) \n",
    "\n",
    "for i, e in enumerate(element):\n",
    "    x1,y1 = node[e][0]\n",
    "    x2,y2 = node[e][1]\n",
    "    x3,y3 = node[e][2]\n",
    "    \n",
    "    A_e = 1/2*np.linalg.det(np.c_[node[e], np.ones((3,1))])\n",
    "    \n",
    "    Phi_e = 1/(2*A_e)*np.array([[x2-x3,x3-x1,x1-x2],\n",
    "                                [y2-y3,y3-y1,y1-y2]])\n",
    "    \n",
    "    if i in Eh:\n",
    "        L1 = np.sqrt((x2-x1)**2+(y2-y1)**2)\n",
    "        L2 = np.sqrt((x3-x2)**2+(y3-y2)**2)\n",
    "\n",
    "        left = np.linalg.inv([[(y1-y2)/L1,(x2-x1)/L1],\n",
    "                              [(y2-y3)/L2,(x3-x2)/L2]])\n",
    "\n",
    "        Phi_L = np.c_[left, np.zeros((2,1))]\n",
    "        \n",
    "        print(Phi_e @ s[e])\n",
    "        Phi_e = np.c_[Phi_e, Phi_L @ Eh[i]]\n",
    "        \n",
    "        ey = np.append(e, -1)\n",
    "\n",
    "        q_e = Phi_e @ s[ey]\n",
    "        print(Phi_e @ s[ey])\n",
    "\n",
    "    else:\n",
    "        q_e = Phi_e @ s[e]\n",
    "    # q_e = Phi_e @ s[e]    \n",
    "    q_h[i] = q_e\n",
    "\n",
    "cfv.figure()\n",
    "\n",
    "# Draw the mesh.\n",
    "\n",
    "cfv.drawElementValues(\n",
    "    ev=q_h[:,1],\n",
    "    coords=coords,\n",
    "    edof=edof,\n",
    "    dofs_per_node=mesh.dofsPerNode,\n",
    "    el_type=mesh.elType,\n",
    "    title=\"Example 01\"\n",
    "        )\n",
    "cfv.showAndWait()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d7b83daa-f7c6-4f98-aeb6-0256c960d087",
   "metadata": {},
   "outputs": [],
   "source": [
    "fp1 = lambda x: -np.cosh(0.4)/np.sin(1)/np.sinh(1)*np.cos(x)\n",
    "fp2 = lambda y: -np.cos(0.6)/np.sin(1)/np.sinh(1)*np.cosh(y)\n",
    "fp3 = lambda x: np.cosh(0.6)/np.sin(1)/np.sinh(1)*np.cos(x)\n",
    "fp4 = lambda y: np.cos(0.4)/np.sin(1)/np.sinh(1)*np.cosh(y)\n",
    "\n",
    "\n",
    "print((fp1(0.6)-fp1(0.4))+(fp2(0.6)-fp2(0.4))+(fp3(0.4)-fp3(0.6))+(fp4(0.4)-fp4(0.6)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d717b115-2bdb-4587-a873-e504a0a5460d",
   "metadata": {},
   "outputs": [],
   "source": [
    "fp1 = lambda x: np.cosh(0)/np.sin(1)/np.sinh(1)*np.cos(x)\n",
    "fp2 = lambda y: np.cos(1)/np.sin(1)/np.sinh(1)*np.cosh(y)\n",
    "fp3 = lambda x: -np.cosh(1)/np.sin(1)/np.sinh(1)*np.cos(x)\n",
    "fp4 = lambda y: -np.cos(0)/np.sin(1)/np.sinh(1)*np.cosh(y)\n",
    "\n",
    "\n",
    "print((fp1(1)-fp1(0))+(fp2(1)-fp2(0))+(fp3(1)-fp3(0))+(fp4(1)-fp4(0)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c8d29a66-eea7-4e82-b2d5-a74a26af9427",
   "metadata": {},
   "outputs": [],
   "source": [
    "from fem import Poisson2D\n",
    "\n",
    "def create_geometry(el_on_curve=20):\n",
    "    # externel boundary ID\n",
    "    dir_id_1 = 100\n",
    "    dir_id_2 = 200\n",
    "    dir_id_3 = 300\n",
    "    neu_id_1 = 400\n",
    "    # hole boundaryd ID\n",
    "    hole_neu_id_1 = 500\n",
    "    hole_neu_id_2 = 600\n",
    "    hole_neu_id_3 = 700\n",
    "    hole_neu_id_4 = 800\n",
    "\n",
    "    # Geometry definition\n",
    "    g = cfg.Geometry()\n",
    "    g.point([0.0, 0.0], ID=0)\n",
    "    g.point([1.0, 0.0], ID=1)\n",
    "    g.point([1.0, 1.0], ID=2)\n",
    "    g.point([0.0, 1.0], ID=3)\n",
    "    g.point([0.4, 0.4], ID=4)\n",
    "    g.point([0.6, 0.4], ID=5)\n",
    "    g.point([0.6, 0.6], ID=6)\n",
    "    g.point([0.4, 0.6], ID=7)\n",
    "    g.spline([0, 1], el_on_curve=el_on_curve, marker=100, ID=0)\n",
    "    g.spline([1, 2], el_on_curve=el_on_curve, marker=101, ID=1) \n",
    "    g.spline([2, 3], el_on_curve=el_on_curve, marker=200, ID=2) \n",
    "    g.spline([3, 0], el_on_curve=el_on_curve, marker=300, ID=3)\n",
    "    g.spline([4, 5], el_on_curve=el_on_curve/4, marker=201, ID=4)\n",
    "    g.spline([5, 6], el_on_curve=el_on_curve/4, marker=301, ID=5)\n",
    "    g.spline([6, 7], el_on_curve=el_on_curve/4, marker=401, ID=6)\n",
    "    g.spline([7, 4], el_on_curve=el_on_curve/4, marker=501, ID=7)\n",
    "    g.surface([0, 1, 2, 3],[[4,5,6,7]])\n",
    "    cfv.drawGeometry(g)\n",
    "    cfv.showAndWait()\n",
    "\n",
    "\n",
    "    bcs = {100:[0],\n",
    "           200:[2],\n",
    "           300:[3],\n",
    "           101:[1],\n",
    "           201:[4],\n",
    "           301:[5],\n",
    "           401:[6],\n",
    "           501:[7],\n",
    "           }\n",
    "\n",
    "    for marker in bcs:\n",
    "        for i in bcs[marker]:\n",
    "            g.curve_marker(ID=i, marker=marker)\n",
    "    cfv.drawGeometry(g)\n",
    "    cfv.showAndWait()\n",
    "    return g\n",
    "\n",
    "\n",
    "def dir100(x,y):\n",
    "    return np.zeros_like(x)\n",
    "\n",
    "def dir200(x,y):\n",
    "    return np.sin(x)/np.sin(1)\n",
    "\n",
    "def dir300(x,y):\n",
    "    return np.zeros_like(x)\n",
    "\n",
    "def neu101(x,y):\n",
    "    return (1/np.tan(1))*np.sinh(y)/np.sinh(1)\n",
    "\n",
    "def neu201(x,y):\n",
    "    return (np.sin(x)/np.sin(1))*(np.cosh(0.4)/np.sinh(1))\n",
    "\n",
    "def neu301(x,y):\n",
    "    return -(np.cos(0.6)/np.sin(1))*(np.sinh(y)/np.sinh(1))\n",
    "\n",
    "def neu401(x,y):\n",
    "    return -(np.sin(x)/np.sin(1))*(np.cosh(0.6)/np.sinh(1))\n",
    "\n",
    "def neu501(x,y):\n",
    "    return (np.cos(0.4)/np.sin(1))*(np.sinh(y)/np.sinh(1))\n",
    "\n",
    "\n",
    "geometry = create_geometry(el_on_curve=20)\n",
    "Dirichlet = [100,200,300]\n",
    "Neumann = [101,201,301,401,501]\n",
    "hole = [[201,301,401,501]]\n",
    "\n",
    "BCfunc = {100:dir100, 200:dir200, 300:dir300,\n",
    "          101:neu101, 201:neu201, 301:neu301, 401:neu401,501:neu501}\n",
    "\n",
    "\n",
    "model = Poisson2D(geometry, Dirichlet, Neumann, hole, BCfunc)\n",
    "model.mesh(showMesh=True)\n",
    "model.efem()\n",
    "print(model.FU)\n",
    "print(np.allclose(model.FU,FU))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a548a8f7-3846-4c8a-8f3f-d7762ebeb581",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "EFEM",
   "language": "python",
   "name": "efem"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
